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Nanoscale materials display enhanced strength and toughness but also larger fluctuations and 
more pronounced size effects with respect to their macroscopic counterparts. Here we study the sys¬ 
tem size-dependence of the failure strength distribution of a monolayer graphene sheet with a small 
concentration of vacancies by molecular dynamics simulations. We simulate sheets of varying size 
encompassing more than three decades and systematically study their deformation as a function of 
disorder, temperature and loading rate. We generalize the weakest-link theory of fracture size effects 
to rate and temperature dependent failure and find quantitative agreement with the simulations. 

Our numerical and theoretical results explain the crossover of the fracture strength distribution be¬ 
tween a thermal and rate-dependent regime and a disorder-dominated regime described by extreme 
value theory. 

PACS numbers: 81.05.ue, 62.20.mm, 62.20.mt, 62.25.Mn 


Nanomaterials have remarkable mechanical properties, 
such as enhanced strength and toughness [1, 2], but dis¬ 
play considerable size effects and sample-to-sample fluc¬ 
tuations, which represent an issue for engineering ap¬ 
plications. Our current understanding of fracture size- 
effects in macroscopic disordered media relies on extreme 
value theory which relates the strength to the statistics 
of the weakest region in the sample [3, 4]. While the 
theory does not consider the effect of stress concentra¬ 
tions and crack interactions, numerical models for the 
failure of elastic networks with disorder show that an ex¬ 
treme value distribution describes failure at large enough 
scales, although the form usually deviates from the stan¬ 
dard Weibull distribution [5-8]. Understanding size ef¬ 
fects in nanomaterials is still an intriguing open issue 
also because of the presence of rate-dependent thermal 
effects that would invalidate the weakest-link hypothe¬ 
sis [9]. Yet, the Weibull distribution is commonly used 
to fit experimental data in carbon based nanomaterials 
[10], although the tensile strength is observed to depend 
on the strain rate [11]. 

Testing fracture properties of graphene is quite chal¬ 
lenging due to the difficulty in applying high tensile 
stresses in a controlled fashion on nanoscale objects [12- 
14]. Therefore numerical simulations represents a vi¬ 
able alternative to understand the size dependence of 
its mechanical behavior [15-19]. Numerical simulations 
of defected carbon nanutubes suggest that failure is de- 
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scribed by the Weibull distribution in quasistatic, zero- 
temperature conditions [20]. Finite temperature molecu¬ 
lar dynamics simulations reveal, however, that the av¬ 
erage tensile strength of nanotubes [21] and graphene 
[19, 22, 23] depends on temperature and loading rate. 
Despite these insightful results, a comprehensive theory 
describing the size dependent fracture strength distribu¬ 
tion of carbon nanomaterials, elucidating the role of ther¬ 
mal fluctuations and strain rate, is still lacking. 

Here we perform large scale molecular dynamics simu¬ 
lations of the deformation and failure of defected mono- 
layer graphene sheets for a wide range of sample sizes, 
vacancy concentration, temperature and strain rate. To 
explain the observed temperature and rate dependence of 
the tensile strength distribution, we generalize extreme 
value theory to the case of thermally activated rate de¬ 
pendent fracture. The resulting theory is shown to be 
in excellent agreement with our simulations and provides 
a general framework to explain rate-dependent thermal 
effects in the failure of disordered nanomaterials. Based 
on our theory, we derive a simple criterion that allows to 
assess the relative importance of structural disorder and 
thermal fluctuations in determining failure. Using this 
rule, one can readily show that the failure of nanoscale 
samples is more prone to thermal induced failure, while 
the fracture macroscopic samples are more likely to be 
ruled by quenched disorder. This confirms previous re¬ 
sults showing that in the limit of very large samples fail¬ 
ure is ruled by extreme value statistics (although not nec¬ 
essarily by the Weibull law) [8, 24]. 

The paper is organized as follows. In section I we de- 
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scribe the molecular dynamics simulation model and in 
section II discuss the numerical results. The theory is 
described in details in section III where we also compare 
its prediction with experiments. Section IV discusses the 
general implications of our work to understand size effects 
in materials at different scales. Appendix A provides de¬ 
tails on the choice of interatomic potential and appendix 
B discusses the fitting method. 


I. MODEL 

We perform numerical simulations of the deformation 
and failure of defected monolayer graphene using the 
LAMMPS molecular dynamics simulator package [25]. 
The carbon-carbon atom interaction is modeled with 
the “Adaptive Intermolecular REactive Bond Order” 
(AIREBO) potential [26]. In order to simulate a realistic 
bond failure behavior, the shortest-scale adaptive cutoff 
of the AIREBO potential has to be fine-tuned [15, 22], as 
detailed in appendix A. The simulated system consists of 
single layer, monocrystalline graphene sheets, composed 
of a variable number N of atoms: N varies from approx¬ 
imately 10^ to 50 X 10^ atoms. The sheets are prepared 
by placing the atoms on a hexagonal lattice; the char¬ 
acteristic lattice length scale A = 1.42 A is chosen so 
that the system is initially in an equilibrium configura¬ 
tion. The sheets have an almost square shape lying on 
the XY coordinate plane; their lateral size depends on N 
and varies between 50 and 360 A (5 and 36 nm). When 
placing defects on the sheets, a fixed fraction of atoms 
is randomly removed; this corresponds to vacancy con¬ 
centrations P = 0.1, 0.2 and 0.5%. While the graphene 
layer is essentially 2D, the atom positions are integrated 
in all the three spatial directions; also, the layers have no 
periodic boundary conditions. 

The simulations are performed by stretching the sam¬ 
ples along the X coordinate axis, corresponding to the 
“armchair” direction of the graphene hexagonal struc¬ 
ture. We select two boundary strips of atoms at the op¬ 
posite X-ends of the sheet. These strips are 3.5 A wide, 
corresponding to 4 atom layers. Hence, the atoms are 
free to move in the Y and Z directions, but follow an 
imposed motion along the stretching direction (X). This 
constraint induces an initial pre-stress on the sheet that is 
visible in the stress-strain curve (see Eig.lb). The Y-end 
boundaries are left free. The system is thermostated by 
means of a Berendsen [27] thermostat with a temperature 
ranging from IK to BOOK, and a characteristic relaxation 
time equal to 0.1 ps; the simulation timestep is set to 
0.5 fs to insure a correct time integration of the atoms 
dynamics. These parameters lead to a slightly under¬ 
damped atom dynamics. Before the stretching protocol 
is started, the system is allowed to relax to thermal equi¬ 
librium from the initial constrained state. Afterwards, 
one of the lateral strips is set in motion, so that the sam¬ 
ple is subject to a constant engineering strain rate 5 inde¬ 
pendent of the system size. The strain rates lie between 
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FIG. 1. Failure of graphene sheets. The graphene sheet 
is composed of A = 50 x 10^ atoms, with a vacancy concen¬ 
tration (porosity) P = 0.1%. The color bar indicates the axx 
component of stress tensor per-atom. a) Graphical view of the 
failure process (from left to right). The crack nucleates from 
one of the defects already present in the material (not neces¬ 
sarily the most stressed) and rapidly grows untill the graphene 
sheet complete failure is achieved, b) The stress strain curve 
displays temperature dependent fracture strength. The pre¬ 
stressed initial condition (e = 0) is due to the constraint ap¬ 
plied to the atoms belonging to the 4 outmost layers of the 
sheet, which are subject to the stretching along X. 


1.28 X 10^s“^ and 1.28 x 10^s“^. As for other molecu¬ 
lar dynamics simulations, the simulated strain rates are 
much higher than those applied experimentally, but the 
deformation speed is still much lower than the sound 
speed in graphene. The chosen strain rate is reached 
by adiabatically ramping up 5, in order to minimize the 
creation of shock waves in the material. As a matter of 
fact, visual inspection of velocity fields shows that shock 
waves are rapidly damped and do not significantly in¬ 
fluence the system dynamics. Simulations are carried 
on until the graphene sheet fractures. Eailure statistics 
are sampled over 100 realizations for each condition in 
which we vary vacancy concentration P, temperature T, 
strain rate 5 and system size N. The only the excep¬ 
tion is provided by systems characterized by T = 300K, 
5 = 0.128 X 10^5-\ A = 20 X 10^ and A = 50 x 10^ 
atoms, where 50 samples were simulated. 


II. SIMULATIONS 

An example of the fracture process is shown in Eig. la, 
where the graphene structure is seen from above at four 
different times during the nucleation and ensuing growth 
of the crack (see also Video 1). The color code represents 
the XX component of the symmetric per-atom stress ten¬ 
sor including both potential and kinetic terms. Typ¬ 
ical stress strain curves are reported in Eig. lb, showing 
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Video 1. The deformation and fracture of a graphene sheet as 
the strain is ramped is shown in the top left panel (P = 0.1%, 
iV = 50 X 10^ T = 300K and e = 0.128 x 10®s“^). The color 
represents the tensile stress <7xx magnitude. A magnihcation 
of the region where the crack is nucleated is shown in the 
bottom left panel. The top right panel reports the same sheet 
viewed under a different angle with a color code representing 
the Z component of the particle positions. The bottom right 
panel reports the corresponding stress strain curve. 


that the tensile strength depends on temperature T. Our 
results provide a clear indication that it also depends on 
system size vacancy concentration P and strain rate 
5, as we discuss below. 

Fig. 2a reports the average failure stress (a) as a function 
of system size for different values of the porosity P, show¬ 
ing that the larger and more defective a sample is, the 
weaker it is. A more complete description of the failure 
statistics is obtained by the survival distribution 5'(cr), 
defined as the probability that a sample has not yet failed 
at a stress a. The numerical results for S{a) are reported 
in Fig. 2b. If a system of volume V fails according to 
extreme value statistics, the survival distribution should 
depend on the volume as S{a) = , where So{(t) 

is the survival distribution of a representative element of 
volume Vb: the smallest independent unit in the sample 
[28]. If we express the volume in terms of the number 
of atoms N and their atomic volume the survival 
probability can be written as S{a) = exp[—A/'Va/Vo/(cr)], 
where f{x) is a suitable function which is a power law x'^ 
in case of Weibull distribution [4], and exponential 
for Gumbel distribution [3]. Fig. 2 shows that the N- 
dependence of the survival distribution follow the pre¬ 
scriptions of extreme value theory, but f{x) is not a 
power law, indicating that the Weibull distribution does 
not represent the data. This is confirmed by the size 
scaling of the average failure stress that does not fol¬ 
low a power law, as would be expected from the Weibull 
distribution. The survival distribution depends also on 
temperature and strain rate, as shown in Fig. 3, which 
is hard to reconcile with the weakest link hypothesis un- 
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FIG. 2. Graphene fracture size effects, a) The average 
failure stress for defected graphene depends on the system 
N size and on the vacancy concentration P. Simulations are 
carried out with T = 300K and £ = 0.128 x 10®The 
lines are the theoretical prediction as discussed in the sup¬ 
porting information. They do not arise as direct fit of the 
numerical curves, but result from the analytical evaluation of 
the integral expression of {a)n- b) The failure stress survival 
distribution at T = 300K,and e = 0.128 x 10®for different 
system sizes with vacancy concentration equal to P = 0.1% 
(blue) , P = 0.2% (green) and P = 0.5% (red). When the 
survival probability distributions are rescaled by N accord¬ 
ing to the predictions of the extreme value theory, the data 
collapse into a single curve that only depends on the vacancy 
concentration P. 


derlying the Weibull distribution. Indeed, by monitoring 
the local stress field axx before failure, we estimate that 
only in less than 20% of the samples (for N = bO x 10^) 
the final crack nucleates in the most stressed region. In 
50-60% of the cases, the final crack is nucleated in re¬ 
gions that ranked fourth of more in terms of stress. This 
is a clear indication that failure is not dictated by the 
weakest link. 
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FIG. 3. Temperature and rate effects of the graphene 
tensile strength distribution. The survival distribution of 
defected graphene sheets (with P = 0.2, £ = 1.28 x 10^ 
and N = 10^) depends on temperature (a) and strain rate (b: 
P = 0.2, T = 300K and N = 10^). The dashed lines represent 
the best least square ht according to the theory of breaking 
kinetics discussed in the text. 


III. THERMAL ACTIVATED FRACTURE OF 
DISORDERED MEDIA 

A. Derivation of the survival distribution 


To understand our simulation results, we generalize 
extreme value theory taking into account thermal fluc¬ 
tuations. We describe the system as a set of n repre¬ 
sentative elements of volume Vb (slabs) such that the 
thermally activated failure of a single element induces 
global failure. Each representative i element obeys lin¬ 
ear elasticity up to a critical strain 5*, so that the elastic 
energy of the sample under an external stress a is given 
by U{cr) = Y.i{Uolei,sl) - Vocrei), where 


Uo{e, 




—oo 


e < Sc 
e > ec, 


( 1 ) 


where E is the Young modulus. The sample is loaded 
at constant strain rate 5 so that a{t) = Eet and critical 
strains are distributed according to a probability density 


function p{Sc)- Assuming the slabs noninteracting and 
identicals, the survival probability for the entire sample 
is given by the product of the survival probabilities of 
each representative element SnicrlT^s) = [So{a\T,s)]^, 
according to the theory of breaking kinetics [29]. The 
representative volume survival probability is defined as 

POO 

S'o((t|T,£) = / dScPiec)'^o{(T\ec,T,e), (2) 

JalE 

where Eq (crj^c, T, 5 ) represents the survival probability 
of a single slab characterized by a failure strain £c- Eq. 
2 reduces to the standard extreme value theory when 
Eo (cr|5c, T, 5 ) = 1, but otherwise depends on tempera¬ 
ture and strain rate. In general, however, the theory pre¬ 
dicts that log(S'^)/n should not depend on the system 
size, as verified by our simulations (see Fig. 2b). 

To estimate the survival distribution of the single slab 
Eq {(j\ec,T,e) we make the phenomenological hypothesis 
that the material failure arises as a thermally activated 
process. Historically, the idea that the solid failure can 
be described by means of the Kramer’s theory, where the 
intrinsic energy barrier is reduced proportionally to the 
applied field, has firstly appeared in material science to 
treat the kinetic fracture of solids under applied stresses, 
and dates back to the works of Tobolsky and Eyring [30] 
and, later, of Zhurkov [31]. More recently it has been 
successfully applied to the study the failure of fibers [32], 
gels [33], wood and fiber-glasses [34] where the potential 
energy barrier is given by the Griffith crack nucleation en¬ 
ergy [35]. Most of previous work focused on the thermal 
dependence of the average strength or the failure time 
in creep experiments and did not address the survival 
distribution and its size dependence. To this end, we 
start from recent theories developed for single-molecule 
pulling, where the molecule rate coefficient for rupture 
(or unbinding) is modified by the presence of an external 
time-dependent force [36-43]. 

In our case, the stress-dependent failure rate of a single 
element characterized by a failure strain 5c, is given by 
an Arrhenius like form [39, 40, 43] 


kia\T,ec) = ko2^^‘^ (^1 - ^"cT] 


( 3 ) 

where ko is the Kramer’s escape rate from the potential 
well described in Eq. 1 [37, 44], 


ko = LUo 


EVo 

ksT 


3/2 



VqEsI 
2kBT ^ 


( 4 ) 


with a characteristic frequency cjo- In onr numerical sim¬ 
ulations one end of the graphene sheet is held fixed, while 
the other is pulled at constant strain rate 5: this can 
be interpreted as the action of a stiff device [40, 43] for 
which Eq.3 has been derived. Eq (cr|5c, T, e) obeys to the 
following first-order rate equation [41] 

rfSo(a|£c,T,£) ^ _kia{t)\T,s,)j:o{a\Sc,T,i), (5) 
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where a{t) = Eet. The survival probability is then read¬ 
ily obtained as 


_ ^ 

Eo {a\sc,T,s) = e 


-(ec-a/E) 

e 


2 VqE 
ksT 


— e 


2 

c k^T 


(6) 

Notice that Eq. 6, only holds for a < ESc since otherwise 
the element fails with probability one (when a Esc 
the Kramer’s theory incorrectly predicts /c(cr|T, a, 5c) — 
0, since it only holds for energy barriers ^ ksT [39]). 
Finally, inserting Eq. 6 in Eq. 2 and, in turns, into the 
constitutive equation for the theory of breaking kinetics, 
we obtain 


Sn{cT\T., 


\Ja/E 


decp{ec)x (7) 


LUO / VqE 


-(s 


— e 


B. Limiting behavior of the theoretical survival 
distribution 


The survival distribution reported in Eq. 7 is written 
as a convolution of the disorder distribution p{sc) with a 
temperature and rate dependent kernel. It is instructive 
to study its limiting behaviors since this allows to assess 
the relevance of thermal and rate dependent effects for 
fracture statistics. Our starting point is the expression 
for the conditional survival probability Eq (cr|5c,T, 5 ) re¬ 
ported in Eq. 6. It is convenient to study its behavior in 
term of the dimensionless parameter A = {V^E)/(kBT)^ 
the ratio between the elastic energy of a representative 
volume element and the thermal energy. In terms of A 
we can write Eo(A) = exp(— G(A)), where 


G(A) 


e 


g-(A£e)"(l-<T/(£c£;))" _ 


(8) 


Thermal fluctuations can be neglected when G(A) —> 0, 
yielding the usual disorder-induced survival probability 
distribution 


Sn{cy\T, 



dScPiSc) 


n 


( 9 ) 


It is interesting to consider first the limit of A ^ oo, 
corresponding to very low temperature and large repre¬ 
sentative volume elements. In this limit, the exponential 
factors in G(A) dominates and the function goes to zero 
even for small strain rates. In more generality, thermal 
fluctuations become negligible when 


5 > UJo^/X 


'^-{Xscfil-a/{seE)f _g-(A£e)" 


( 10 ) 


Another interesting limit is the low stress regime (i.e. 
^ > 0) where 


^0 (cp|5c, T, 5 ) ^ 1 


2^ ^ 


3/2 


5cCr. 


( 11 ) 


Hence, thanks to Eq. 2, the survival distribution for a 
representative element is given by 


So (c^|T, 5 ) ^ 1 


^ [E ( 

e V Fo \kBT) 


{£c)o- 


( 12 ) 


where (sc) = decScp(sc)- Therefore, the survival 
probability distribution function for the entire system can 
be recast as 


lnSn{< 7 \T,s) ^LUO [E ( ^ 

- n - 


(13) 


displaying a linear dependence on the applied stress, irre¬ 
spective of the failure strain distribution function p(5c). 


C. Fit of the numerical data 

Eq. 7 provides an excellent fit to the results obtained 
from numerical simulations of defected graphene at dif¬ 
ferent defect concentrations P, temperature T and load¬ 
ing rate 5 . To fit the numerical simulations with Eq. 
9, we first need to establish the form of p(5c). This 
is a phenomenological function describing the distribu¬ 
tion of failure strains of representative volume elements 
at zero temperature. A reasonable estimate of its func¬ 
tional form can be obtained from simulations at low tem¬ 
perature (i.e. T = IK), where thermal fluctuations are 
negligible, as discussed in details in appendix B. The nu¬ 
merical outcomes indicate that p{sc) follow the Gumbel 
distribution [3] (see Fig. 4). We then insert the resulting 
form of p{ec) in Eq. 7 which we adopt as a fitting func¬ 
tion for the numerical survival probability 5'(cr), with ujq 
and Vq as fitting parameters. 

The representative volume Vq ranges between O.lnm^ 
and 0.25nm^, while the characteristic frequency is found 
in the range ^ 6 x 10^5“^ and ^ 10^5“^ (see Fig. 5) 
Moreover, from the survival distribution we also calcu¬ 
late, without additional fitting, the system size depen¬ 
dence of the average tensile strength (cr)^, which displays 
an excellent agreement with simulations results as shown 
in Fig. 2a. Further details on the fitting methodology 
and the analytical expressions used in our model are re¬ 
ported in appendix B. 


IV. DISCUSSION 


Therefore there is a temperature and stress dependent 
critical strain rate above which we can neglect thermal 
fluctuations. 


In conclusions, we have performed extensive numerical 
simulations for the tensile failure of defected graphene fo¬ 
cusing on the size effects of the strength distribution for 
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FIG. 4. Survival distribution of defected graphene at 
low temperature. We report the survival distribution ob¬ 
tained from simulations at T = IK, £ = 0.128 x 10® 
and A/" = 10^ for different values of the vacancy concen¬ 
tration P. The numerical data are fitted with the expo- 

cr 

nential function Ae^^o (solid lines), leading to a Gumbel 

£c 

£c _£Q 

distribution for the failure strains p{sc) = Ae^o (see 

Eq.Bl). For P = 0.1% we obtain A = 7.92 ± 0.05 x 10“®®, 
£0 = 0.00125 ± 0.00004. For P = 0.2%: A = 1.767 ± 
0.005 X 10“®®, £0 = 0.001338 ± 0.000007. For P = 0.5% 
A = 1.804 ± 0.007 X 10“^®, £o = 0.00167 ± 0.00004. 


FIG. 5. Fitting parameters for graphene survival dis¬ 
tributions. The best fitted values of a) the representative 
volume Vo and b) the activation frequency as a function of 
temperature {N = 10^, £ = 1.28 x 10®s“^). These values 
are obtained by the least square fit of the numerical survival 
probability distribution with the expression B2 (dashed lines 
in Fig.3a). The same values, c) and d), as a function of the 
strain rate {N = 10^, T = 300K), obtained from the best fits 
shown in Fig.3b (P = 0.2%) and Fig.7 (P = 0.5%) (dashed 
lines). 


different temperatures and loading rate. The results of 
numerical simulations show deviations from the weakest- 
link hypothesis but can be explained by taking explicitly 
into account the effect of thermally activated crack nu- 
cleation. The resulting theory describes well our results 
and could prove useful to understand the tensile strength 
distribution of other nanomaterials such as carbon nan¬ 
otubes or other nanowires. 

At present it is not possible to compare our numeri¬ 
cal and theoretical predictions directly to experiments. 
Experimental measurements of the strength of graphene 
sheets are mostly based on indentation tests [12], while 
tensile tests only recently appeared in the literature [14] 
but thermal, rate and size effects have not been stud¬ 
ied. Furthermore, most experimental studies are focus¬ 
ing on the strength of graphene in pristine conditions 
[12] without defects or pre-existing cracks. Our theory 
is, however, very general yielding predictions that should 
be applicable also to other carbon nanomaterials and al¬ 
lows to formulate general considerations on the relevance 
of thermal effects for fracture. 

Eq. 7 suggests that thermal fluctuations can be ne¬ 
glected for large enough strain rate, since in this limit 
Eo — 1 and the sample fails according to the weakest 
link statistics. In our simulations, we have E ^ lO^^Pa, 
Vo — O.lnm^ so that at room temperature we estimate 
A 10®. If we use this value in Eq. 8, we find that 
for 5c — 0.1 the exponential terms do not vanish close 


to failure (i. e. for a > 0.9{Esc)) and thermal effects 
should therefore be relevant. Indeed using ^ 10“^ in 

Eq. 10, one can readily show that thermal effects start 
to become relevant for T > lOK in agreement with our 
simulations. 

The same argument suggests that in macroscopic sam¬ 
ples, with larger representative volume elements, ther¬ 
mally activated failure can often be ignored, even at room 
temperature. Consider for instance a ceramic material, 
like sintered a-alumina [45], with E = lO^^Pa and a typi¬ 
cal tensile strength of cr = 10^Pa. Assuming that the rep¬ 
resentative volume element corresponds to a grain size of 
Vo — l(/im)^, we can estimate A 10^^. Now the expo¬ 
nential factors impose that G(A) ^ 0 even at low strain 
rates, implying that the strength distribution should be 
described by conventional extreme value theory. Indeed, 
experiments show that the strength distribution is de¬ 
scribed by Weibull statistics with parameters that are 
largely temperature independent [45]. Our theory thus 
provides a simple way to estimate the relevance of ther¬ 
mal and rate dependent effects for fracture. This re¬ 
sult could have important implications for applications 
to micro- and nano-mechanical devices whose reliability 
may crucially depend on the control of thermally acti¬ 
vated failure. 
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Appendix A: Interatomic potential and cutoff tuning 


The carbon-carbon atom interactions were modeled 
using the “Adaptive Intermolecular REactive Bond Or¬ 
der” (AIREBO) potential [26], which was originally de¬ 
veloped as an extension of the “REactive Bond Order” 
potential (REBO) [46]. In turn, the REBO potential 
was developed to describe covalent bond breaking and 
forming with associated changes in atomic hybridization 
within a classical potential; it has proven an useful tool 
for modelling complex chemistry in large many-atom sys¬ 
tems. The AIREBO potential improves the REBO po¬ 
tential with an adaptive treatment of non-bonded and 
dihedral angle interactions that is employed to capture 
the bond breaking and bond reformation between carbon 
atom chains. The analytical form of the AIREBO poten¬ 
tial (as discussed in the documentation [25]) is written 
as: 


BEBO 

The term has the same functional form as 

the hydrocarbon REBO potential developed in [46]. We 
will not cover here the details of the energetic terms 
which are thoroughly discussed in the mentioned ref¬ 
erence. In short, the REBO term gives the model its 
short to medium range reactive capabilities, describing 
short-ranged C-C, C-H and H-H interactions (r < 2 A). 
These interactions have strong coordination-dependence 
through a bond order parameter, which adjusts the at¬ 
traction between the i^j atoms based on the position 
of other nearby atoms and thus has 3- and 4-body de¬ 
pendencies. A more detailed discussion of formulas for 
this part of the potential are given in [26]. The 
term adds longer-ranged interactions (2 < r < ^^^toff 
A) using a form similar to the standard Lennard-Jones 
potential. It contains a series of switching functions so 
that the short-ranged LJ repulsion (1/r^^) does not inter- 
fere with the energetics captured by the term. 

The extent of the e\^^ interactions is determined by a 
cutoff argument; in general the resulting E\j^ cutoff is 
approximately 10 A, in this work we consider a cutoff 
of approximately 14 A. Einally, the term is an 

explicit 4-body potential that describes various dihedral 
angle preferences in hydrocarbon configurations. 

The AIREBO potential has been extensively used to 
simulate and predict mechanical properties of carbon- 
based materials, i.e. fullerene, carbon nanotube and 
graphene [22]. Eurthermore, it offers a valid tradeoff be¬ 
tween accuracy and computational efficiency; a realistic 
fracture of large system sizes can be simulated in reason¬ 
ably short time scales (a few hours on recent computers). 
Other interaction models can offer little improvement to 


the actual realism of the simulation, at the cost of much 
larger computational costs: for example, the ReaxEE po¬ 
tential, or DET semiclassical approaches could describe 
more accurately the fast time scales of chemical reactions, 
but this would not change the ultimate failure length of 
the C-C bond: the expected maximum elongation for a 
C-C bond in graphene is around 0.178 nm. On the other 
hand, the use of faster but too simplistic models (e.g. 
Lennard-Jones potentials, mass and spring systems, or 
other elastic models) fail to significantly reproduce a re¬ 
alistic behavior. 

However, in order to simulate a realistic bond failure 
behavior, the short-scale C-C adaptive cutoff (rc) of the 
AIREBO potential has to be tuned. In fact, it has been 
observed [15, 47] that, during simulations of fracture of 
covalent bonds and without cutoff tuning, the shortest- 
scale potential introduces a sharp increase of bond forces 
near the cutoff distances, which in turn causes spurious 
increase in fracture stress and strain [22]. It should also 
be noted here that this phenomenon is specifically rele¬ 
vant for perfect graphene and CNT lattices, while it is 
much less pronounced in defected samples, due to the 
disorder induced in the lattice by the atom vacancies. 
This issue has been solved in the past by incrementing 
the short-scale cutoff lenght of the potential; the cited 
papers increase this parameter to 2.0 A. This, however, 
has the side effect of leading to a singular behavior in 
the atomic pair potential when the atom atom distance 
is exactly 2 A. 



FIG. 6. Tuning the AIREBO potential. The stress-strain 
curve obtained as a function of the cutoff rc. The simulated 
graphene sheet is composed by A = 10^ carbon atoms at 
T = 300K, without in-built defects (P = 0) and pulled at 
constant strain rate e = 0.128 x 10®For Tc < 0.195nm, 
the stress displays a spurious increase while for Vc = 2nm the 
pair potential shows an unphysical singularity (not shown). 
The chosen value of Vc is set to 0.195nm. 

We performed stretching simulations varying the cut¬ 
off parameter from Vc = 0.17 nm (default value) to 
Vc = 0.2 nm in both armchair (X) and zigzag (Y) direc- 










tions of the graphene sheet with no vacancies (P = 0). 
The stress-strain curve obtained from the numerical sim¬ 
ulations are shown in Fig. 6 . For Vc < 0.195 nm, a 
sharp increase on tensile stress for large strains is ob¬ 
served, leading to an unphysical ultra-high failure stress 
and corresponding failure strain. Increasing the Vc in 
the range 0.195 < Vc < 0.2 nm strongly suppresses this 
phenomenon. Moreover, the stress-strain data reported 
in Fig. 6 , clearly display that the failure strain varies 
from 0.13 to 0.25 when Vc is in the range 1.95 < Tc < 2, 
whereas the failure stress exhibites a much weaker fluc¬ 
tuation (from 85 x 10^ Pa to 95 x 10^ Pa). Finally we 
notice that for defected samples like those investigated in 
the present article, i.e. P 7 ^ 0, the values of the failure 
stresses and strains do show a much less marked depen¬ 
dence on the choice of whenever 1.95 < Vc < 2. 

Appendix B: Details of the fitting method 

To fit the numerical simulations with Eq. 7, we first 
obtain p{ec) from simulations at low temperature (i.e. 


T = IK). As shown in Fig.4, the numerical survival 
distribution function — , obtained at T = IK and 

5 = 0.128x 10^5-\ can be nicely fitted with the following 
exponential form Ae~^. The theoretical prediction for 
the survival probability distribution furnished by Eq.9 
requires —In dscp{£c) = Ae~^ ^ once we assume 
that Vq = Va when T ^ 0. Hence, we obtain 

p{ec) = , (Bl) 

which corresponds to a Gumbel distribution of failure 
strains [3]. The numerical values of the fitting parame¬ 
ters A, So are reported in the caption of Eig.4 for three 
vacancy concentrations P. We notice that the simulated 
samples for T = IK are 250 in the case of P = 0.1% and 
850 for P = 0.2% and 800 P = 0.5%. 

We then perform the least square fit of the numerical sur¬ 
vival probabilities — , obtained for different values 

of T, 5 and porosities P (see Eig.s 3a,b, 7, 8), with the 
following function 


hiSn{(j\T,s) 

N 


Ki 


dScAe^o 


£c_Ae^o 


VqE 


j/E 


-{ec-cT / E) 


2 VqE 
knT _ 


,2 VqE 
^ knT 


(B2) 


where the fitting parameters are the representative vol¬ 
ume element Vo and the characteristic frequency ujo • The 
atomic volume 14 has been evaluated by considering a 
density of 38.18 atoms per nm? and a sheet thickness 
equal to 0.335 nm, yielding I 4 = 8.744 x 10“^nm^. Eor 
any value of P, the corresponding values of A and So ob¬ 


tained from the best fit of the data in Eig.4 are plugged 
into Eq.B2. The fitted Vo and uJo corresponding to Eig.s 
3a,b, 7, 8 are reported in Eig.5. 

Einally we provide the analytical expression for the 

distribution of failure stresses defined as Pn(cr|T, 5 ) = 

_ dS-^ , 

da 


E 


+ 


Pn (o-|r,e) = n 

^ii^) ' {sc-%)e 


a- —SO. / ^qE 

^ - Ae'^ ^ 

Ae^^o e 


_Vqs 1 

l_e Ek^i 


^2 VqE 

T g 


^0 / VqE 


+ 


,>,2 VqE _2 ^ qE 

ksT _g ^ kjgT 


(B3) 


Eq.B3 allows to derive the mean failure stress as 
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FIG. 7. Survival distribution of defected graphene at 

P=0.5%. We report the survival distribution obtained from 
simulations at T = 300K, iV = 10^ and three strain rates 
e. The vacancy concentration is set to P = 0.5%. The fit¬ 
ting function is provided by Eq.B2: the values of the fitting 
parameters Vb and uq are reported in Fig.5. 


{a)n = n /y da ^Sn-i {cr\T,e) < 


A - -Ac 

Ae^^o e 


a 

Eeo ^ Y ' 


VqE 




l_e Ek^i 


+ 


+ ' pEd^cAe^o (ee-i)e ^ 


_ ^ / ^qE 




3 (^c ^/E)- k^T ‘^ck^T 


(B4) 


> . 


This quantity can be analytically calculated and plotted 
as a function of setting n = as shown in Fig.2a 

for T = 300K, 5 = 0.128 x 10^5“^ and three values of 
the vacancy concentration P. We emphasize that in this 
case no fit, but just the numerical evaluation of the in¬ 
tegral expression of {(j)n B4 is provided, making use of 
the proper values of ^4,5o, Vb, cjo, obtained by fitting the 
survival probabilities displayed in Fig.4 and Fig.8. 
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[GPa] 


FIG. 8. Survival distribution of defected graphene at 
different P We report the survival distribution obtained 
from simulations at T = 300K, £ = 0.128 x 10®and 
A/" = 10^ for three vacancy concentrations P. The numeri¬ 
cal data are fitted with the expression B2 (dashed lines) us¬ 
ing the proper values of A and eo reported in Fig.4. The 
fitted values of Vb and ujq are shown in Fig.5 for P = 0.2% 
and P = 0.5%. For P = 0.1% the least square fit gives 
Vb = 0.3806 ±0.0003nm® and c^o = 4.1416 ± 0.0006 x 10^s“^ 
The set of parameters A,£o,Vo and c^o which characterize 
uniquely the theoretical expression B2 (dashed lines) are in¬ 
serted into Eq.B4 to calculate the mean average rupture stress 
{a)n as a function of N, shown in Fig.2a. 
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